---
title: "Replication Material -- Cutting trees to balance budgets: the IMF's role in deforestation --"
author: "Forster, T., Bhandary, R.R. & Gallagher, K.P."
date: "December 2025"
output:
  pdf_document: default
  html_document: default
---

# Set up

```{r}
# Load packages
library(dplyr)
library(magrittr)
library(tidyr)
library(ggplot2)
library(scales)
library(quanteda)
library(stringr)
library(stargazer)
library(patchwork)
library(hrbrthemes)
library(cowplot)
library(modelsummary)
library(clubSandwich)
library(fixest)
library(purrr)
library(broom)
```


```{r}
# Set path/directory
dir_data <- "~/Replication/"
```




# Text analysis: Conditionality

```{r}
# Load data: INET Conditionality, 1980-2019
# Source: imfmonitor.org
dt_cond_raw <- read.csv(paste0(dir_data, "/IMFMonitor2023_Conditions_Raw.csv"), 
                        header = TRUE) %>% 
  mutate(Condition.ID = row_number())
```


## Text pre-processing

We first clean the texts and convert the data into a text corpus, defining `text` as the text-variable and using `row_id` as the identifier. This includes tokenizing into unigrams. In addition, we use lemmatization to achieve word equivalence. Before these steps, we remove the genitive s. 

```{r}
# Remove "*'s"
dt_cond_edit <- dt_cond_raw %>% 
  mutate(text = str_replace_all(Condition.Text, pattern = "'s", replacement = ""))

# Corpus
corp_cond <- corpus(dt_cond_edit, 
                    text_field = "Condition.Text", 
                    docid_field = "Condition.ID")

# Tokenization
tok_cond <- tokens(corp_cond,
                   remove_punct = TRUE, remove_symbols = TRUE, 
                   remove_separators = TRUE, split_hyphens = TRUE,
                   include_docvars = TRUE)

# Lemmatization
list_lemma <- read.csv(paste0(dir_data, "/lemmatization-en.txt"), 
                       sep = "\t", as.is = TRUE, 
                       header = FALSE,
                       fileEncoding = "UTF-8-BOM")
lemma_cond <- tokens_replace(tok_cond, 
                             pattern = list_lemma[, 2], 
                             replacement = list_lemma[, 1],
                             valuetype = "fixed", case_insensitive = TRUE)

# Convert to dfm
dfm_cond <- dfm(lemma_cond)
```


Strictly speaking, we do not need the DFM object because we apply the dictionary (see next steps) to the token object.

```{r}
dfm_cond
nfeat(dfm_cond)
ndoc(dfm_cond)
```

Let's save the output.

```{r}
saveRDS(lemma_cond, paste0(dir_data, "/lemma_IMFConditions_Dec25.rds"))
saveRDS(dfm_cond, paste0(dir_data, "/dfm_IMFConditions_Dec25.rds"))
```


## Dictionary

### Define parameters

```{r}
lemma_cond <- readRDS(paste0(dir_data, "/lemma_IMFConditions_Dec25.rds"))
dfm_cond <- readRDS(paste0(dir_data, "/dfm_IMFConditions_Dec25.rds"))
```


Let's check some keywords in context to see what these terms are.

```{r}
# Not run
# kwic(lemma_cond, pattern = c("forest", "tree"))
```


We use the following keywords to approximate forest conditionality.

```{r}
# Define dictionary
forest <- c("forest", "forests", "rainforest", "rainforests", "wood", "woods", "forestry", "logging", "felling", "deforest", "deforestation", "desertify", "desertification", "tree", "trees")

# Dictionary object
dict_forest <- dictionary(list(Forest = forest))
```


### Apply dictionary

Now we apply the dictionary to the token object with `tokens_lookup()`.

```{r}
# Dictionary DFM
tok_dict <- tokens_lookup(lemma_cond, valuetype = "fixed",
                               dictionary = dict_forest,
                               case_insensitive = FALSE)

# Convert tokens object into DFM
dfm_dict <- dfm(tok_dict)

# Reshape data
df_dict_output <- convert(dfm_dict, to = "data.frame") %>% 
  rename(Condition.ID = doc_id) %>%
  mutate(Condition.ID = as.numeric(Condition.ID)) %>% 
  # add meta data
  left_join(dt_cond_edit, by = c("Condition.ID"))

# Save output
saveRDS(df_dict_output, paste0(dir_data, "/dict_output_Dec25.rds"))
```


After manually inspecting all the positives (i.e., conditions identified as pertaining to forest conditionality), we fix one false positive related to Bretton Woods.

```{r}
df_dict_output <- readRDS(paste0(dir_data, "/dict_output_Dec25.rds")) %>% 
  # Manual fix of Bretton Woods
  mutate(forest = if_else(Condition.ID == 8019, 0, forest))
```


```{r}
# Aggregate to Country-Year
dt_FOR <- df_dict_output %>% 
  filter(forest >= 1) %>% 
  group_by(Country.Code, Condition.Year) %>% 
  summarize(BA1FOR = n(), 
            .groups = "keep") %>% 
  ungroup()
```


Attach this to the original dataframe with conditionality --- and aggregate to the types of conditionality as discussed in text.

```{r}
# IMF Conditionality data
# Source: imfmonitor.org
dt_INET_raw <- haven::read_dta(paste0(dir_data, "/IMFMonitor2023_Conditions_Main.dta"))
```


```{r}
# Create IMF condition counts
dt_INET_new <- dt_INET_raw %>% 
  left_join(dt_FOR, by = c("ccode" = "Country.Code",
                           "year" = "Condition.Year")) %>%
  # Group policy areas
  mutate(BA1NEW = BA1INS + BA1PRI + BA1SOE + BA1LAB + BA1SP + BA1POV,
         BA1FISCAL = BA1DEB + BA1RTP + BA1FP) %>% 
  # select vars
  select(year, BA1TOT, BA1ENV, BA1FOR,
         BA1FISCAL, BA1EXT, BA1FIN, BA1NEW)
```



# Figures

## Figure 1


```{r}
visualize_FOR <- dt_INET_raw %>% 
  left_join(dt_FOR, by = c("ccode" = "Country.Code",
                           "year" = "Condition.Year")) %>%
  # new groups
  mutate(BA1NEW = BA1INS + BA1PRI + BA1SOE + BA1LAB + BA1SP + BA1POV,
         BA1FISCAL = BA1DEB + BA1RTP + BA1FP) %>% 
  # select vars
  select(year, BA1TOT, BA1ENV, BA1FOR,
         BA1FISCAL, BA1EXT, BA1FIN, BA1NEW) %>% 
  group_by(year) %>% 
  summarize(across(everything(), .fns = \(x) sum(x, na.rm = TRUE)))
```


```{r}
p <- visualize_FOR %>% 
  ggplot(aes(x = year)) +
  geom_col(aes(y = BA1FOR), fill = "skyblue") +
  geom_point(aes(y = BA1ENV), color = "#444444") +
  geom_line(aes(y = BA1ENV), color = "#444444", linewidth = 0.7) +
  annotate("text", x = 1989, y = 14, 
           label = "Environmental \nreforms", color = "#444444",
           family = "Arial", size = 2.5) +
  annotate("text", x = 2006, y = 8, 
           label = "Forest \nmanagement", color = "skyblue",
           family = "Arial", size = 2.5) +
  hrbrthemes::theme_ipsum_rc(base_family = "Arial") +
  theme(legend.position = "None",
        axis.title.x = element_text(size = 9),
        axis.title.y = element_text(size = 9),
        axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8)) +
  labs(x = "Year",
       y = "Number of conditions") + 
  scale_y_continuous(breaks = c(0, 8, 16, 24),
                     limits = c(0, 25))

# Save at journal size: 1.5 columns = 114 mm wide
ggsave(
  filename = file.path(dir_data, "/Figure1_Dec25.jpg"),
  plot = p,
  width = 114,
  height = 71.25,
  units = "mm",
  dpi = 720,
  device = ragg::agg_jpeg,
  bg = "white"
)
```


## Figure 2

All conditions with groups used below.

```{r}
viz_stacked <- visualize_FOR %>% 
  # drop BA1TOT, BA1FOR
  select(-c(BA1TOT, BA1FOR)) %>% 
  # pivot long
  pivot_longer(cols = c("BA1ENV", "BA1FISCAL", "BA1EXT", "BA1FIN", "BA1NEW"),
               names_to = "Policy.Area",
               values_to = "No.Conditions") %>% 
  # adjust labels
  mutate(Policy.Label = case_when(
    Policy.Area == "BA1ENV" ~ "Environment",
    Policy.Area == "BA1FISCAL" ~ "Fiscal policy",
    Policy.Area == "BA1EXT" ~ "External sector",
    Policy.Area == "BA1FIN" ~ "Monetary policy",
    Policy.Area == "BA1NEW" ~ "Non-core issues")) %>% 
  mutate(Policy.Label = factor(Policy.Label, 
                               levels = c("Environment", "Fiscal policy", "External sector", "Monetary policy", "Non-core issues"),
                               ordered = TRUE))
```


```{r}
p2 <- viz_stacked %>%
  ggplot(aes(x = year, y = No.Conditions, fill = Policy.Label)) +
  geom_bar(stat = "identity") +

  hrbrthemes::theme_ipsum_rc(base_family = "Arial", base_size = 8) +

  theme(
    legend.position = "bottom",
    legend.title.align = 0.5,

    # Axis text & titles (≥ 8 pt)
    axis.title.x = element_text(size = 8, family = "Arial", margin = margin(t = 4)),
    axis.title.y = element_text(size = 8, family = "Arial", margin = margin(r = 4)),
    axis.text.x  = element_text(size = 8, family = "Arial"),
    axis.text.y  = element_text(size = 8, family = "Arial"),

    # Legend text at 6 pt
    legend.title = element_text(size = 6, family = "Arial", margin = margin(b = 2)),
    legend.text  = element_text(size = 6, family = "Arial"),

    # Make legend box tighter
    legend.key.size  = unit(3, "mm"),
    legend.spacing.x = unit(2, "mm"),
    legend.spacing.y = unit(1, "mm"),
    legend.margin    = margin(t = 2, b = 2, l = 2, r = 2),
    legend.box.margin = margin(t = -6),

    # Shrink outer plot margins
    plot.margin = margin(t = 4, r = 4, b = 2, l = 6)
  ) +

  guides(
    fill = guide_legend(
      title.position = "top",
      title.hjust    = 0.5,
      nrow = 2,
      byrow = TRUE
    )
  ) +

  labs(
    x = "Year",
    y = "Number of conditions",
    fill = "Policy area"
  ) +

  scale_y_continuous(labels = scales::comma)

ggsave(
  filename = file.path(dir_data, "/Figure2_Dec25.jpg"),
  plot = p2,
  width = 85,
  height = 53.1,
  units = "mm",
  dpi = 720,
  device = ragg::agg_jpeg,
  bg = "white"
)
```



Some diagnostics. 1,991 observations with IMF programs (`BA1TOT > 0`).

```{r}
dt_INET_new %>% 
  filter(BA1TOT > 0)
  # filter(BA1FISCAL >= 10)
  # filter(BA1FISCAL >= 1)
  # filter(BA1FIN >= 1)
```

1,422 (71.4%) with FISCAL (`filter(BA1FISCAL >= 10)`);
1,956 (98.2%) with EXT `filter(BA1FISCAL >= 1)`;
1,894 (95.1%) with FIN `filter(BA1FIN >= 1)`



# Regression Analysis

## Sample

## Data

```{r}
dt_sample_1pct <- read.csv(paste0(dir_data, "/IMFFOR_sample.csv"),
                      header = TRUE, fileEncoding = "UTF-8-BOM") %>%
  group_by_all() %>% 
  expand(year = c(2000:2019)) %>% 
  ungroup()
```


### Dependent Variable

First, let's look at the dependent variable. Log of net annual (non-fire related) tree cover loss (in ha) from Tyukavina et al. (2022); available on the website of the World Resources Institute.

```{r}
# Forest loss net of fires
dt_tcloss_WRI <- read.csv(paste0(dir_data,
                                     "/treecover_loss_from_fires__ha_GFW_20251214.csv"),
                              fileEncoding = "UTF-8-BOM", header = TRUE) %>%
  filter(Country.Name_WDI != "NA") %>% 
  rename(ccode_iso3 = iso_GFW,
         year = umd_tree_cover_loss__year,
         umd_tcloss_ha = umd_tree_cover_loss__ha,
         umd_tcloss_fire_ha = umd_tree_cover_loss_from_fires__ha) %>% 
  mutate(umd_lnnettcloss_ha = log(umd_nettcloss_ha + 1),
         L1.year = year - 1) %>% 
  select(Country.Name_WDI, L1.year, umd_tcloss_ha, umd_lnnettcloss_ha)
```

```{r}
summary(dt_tcloss_WRI$umd_lnnettcloss_ha)
```

### Treatment variable: IMF conditionality

Now we load the IMF Conditionality data from Kentikelenis & Stubbs (2023), available from imfmonitor.org. We aggregate these condition counts (at the level of the country-year).

```{r}
dt_INET_raw <- haven::read_dta(paste0(dir_data, 
                                      "/IMFMonitor2023_Conditions_Main.dta")) %>% 
  mutate(IMFBA1 = if_else(BA1TOT > 0, 1, 0),
         IMFBA2 = if_else(BA2TOT > 0, 1, 0),
         BA1FISCAL = BA1DEB + BA1RTP + BA1FP,
         BA2FISCAL = BA2DEB + BA2RTP + BA2FP,
         BA1NEW = BA1INS + BA1PRI + BA1SOE + BA1LAB + BA1SP + BA1POV,
         BA2NEW = BA2INS + BA2PRI + BA2SOE + BA2LAB + BA2SP + BA2POV) %>%
  rename(ccode_INET = ccode) %>% 
  # select vars
  select(ccode_INET, year, 
         IMFBA1, IMFBA2,
         BA1TOT, BA1ENV, BA1FISCAL, BA1EXT, BA1FIN, BA1NEW,
         BA2TOT, BA2ENV, BA2FISCAL, BA2EXT, BA2FIN, BA2NEW) %>% 
  # remove NA
  filter(ccode_INET != "NA")
```


Further, we attach a binary indicator for program participation, IMF55 (if an IMF program has been active for at least five calendar months per year; derived from Kentikelenis & Stubbs 2023).

```{r}
# binary indicator
dt_IMF55_raw <- read.csv(paste0(dir_data, "/IMF55_1980-2019.csv"),
                         header = TRUE, fileEncoding = "UTF-8-BOM") %>% 
  select(-ccode) %>% 
  # remove NA
  filter(cname != "NA") 
```



Here, calculate number of years under an IMF program for the V20 countries.

```{r}
vec_V20_AFR <- c(
  # Africa & Middle East
  "Burkina Faso",
  "Benin",
  "Chad",
  "Comoros",
  "Cote d'Ivoire",
  "Congo, Dem. Rep.",
  "Ethiopia",
  "Eswatini",
  "Gambia",
  "Ghana",
  "Kenya",
  "Lebanon",
  "Liberia",
  "Madagascar",
  "Malawi",
  "Morocco",
  "Mozambique",
  "Namibia",
  "Niger",
  "Palestine",
  "Rwanda",
  "Senegal",
  "Sierra Leone",
  "South Sudan",
  "Sudan",
  "Tanzania",
  "Togo",
  "Tunisia",
  "Uganda",
  "Guinea",
  "Jordan",
  "Yemen")

vec_V20_AP <- c(
  # Asia-Pacific
  "Afghanistan",
  "Bangladesh",
  "Bhutan",
  "Cambodia",
  "Fiji",
  "Kiribati",
  "Kyrgyz Republic",
  "Maldives",
  "Marshall Islands",
  "Mongolia",
  "Nepal",
  "Palau",
  "Papua New Guinea",
  "Pakistan",
  "Philippines",
  "Samoa",
  "Sri Lanka",
  "Timor-Leste",
  "Tonga",
  "Tuvalu",
  "Vanuatu",
  "Vietnam")


vec_V20_LAT <- c(
  # Latin America & Caribbean
  "Barbados",
  "Colombia",
  "Costa Rica",
  "Dominica",
  "Dominican Republic",
  "Grenada",
  "Guatemala",
  "Guyana",
  "Haiti",
  "Honduras",
  "Nicaragua",
  "Paraguay",
  "Saint Lucia",
  "Trinidad and Tobago")

vec_V20 <- c(vec_V20_AFR, vec_V20_AP, vec_V20_LAT)
```


```{r}
dt_IMF55_raw %>%
  filter(year >= 1980 & year <= 2019) %>% 
  filter(cname %in% vec_V20) %>% 
  filter(IMF55 == 1) %>% 
  select(cname) %>% 
  unique()
```

53 countries

852 years with an IMF program

53 * 40 = 2120

i.e., 40.2%


### Controls

First, let's load the WDI data --- the World Bank's World Development Indicators (WDI) provide data on key macroeconomic and social indicators. We extracted selected variables from the latest version (July 2025): https://datatopics.worldbank.org/world-development-indicators/

```{r}
dt_WDI_var <- readRDS(paste0(dir_data, "/dt_WDI_Dec2025.rds"))
```


Next, we also collect data on Current account imbalance (% GDP) [BCA_NGDPD] from the IMF WEO, Version October 2025.

This version only records NAs for Syria. Yet historic data (https://www.imf.org/external/datamapper/BCA_NGDPD@WEO/SYR?zoom=SYR&highlight=SYR) does include data for Syria, which we add manually below.

```{r}
# WEO Oct 2025
dt_WEO_raw <- read.csv(paste0(dir_data, "/dataset_T20_23_57.803850913Z_IMF.RES_WEO_9.0.0.csv"),
                   fileEncoding = "UTF-8-BOM") %>% 
  rename(Country.WEO = COUNTRY,
         year = TIME_PERIOD,
         cabal_WEO = OBS_VALUE) %>% 
  select(Country.WEO, year, cabal_WEO)

# Manual file for Syria
dt_WEO_Syria <- tibble(
  Country.WEO = rep("Syrian Arab Republic", n = 11),
  year = c(2000:2010),
  cabal_WEO = c(5.2, 4, 4.9, -7.7, -3.1, -2.2, 1.4, -0.2, -1.3, -2.9, -2.8)
)

# WEO final
dt_WEO <- dt_WEO_raw %>% 
  rbind(dt_WEO_Syria)
```


Third, we extract data from VDem, using the `vdemdata` package.

```{r}
# Vdem controls
dt_vdem_raw <- vdemdata::vdem

# subset, final
dt_vdem <- dt_vdem_raw %>% 
  filter(year >= 2000 & year <= 2025) %>%
  select(country_text_id, year, 
         # VDem
         v2x_libdem, v2xel_elecparl, v2xel_elecpres,
         # WGI
         e_wbgi_pve) %>%
  rename(libdem_VDem = v2x_libdem,
         legelec_VDem = v2xel_elecparl,
         exelec_VDem = v2xel_elecpres,
         polstab_WGI = e_wbgi_pve) %>%
  # Dummy for any elections
  mutate(elec_VDem = pmax(legelec_VDem, exelec_VDem))
```


Fourth, we add data on UNGA voting affinity.

```{r}
# UNGA
dt_unga <- readRDS(paste0(dir_data, "/UNGAvoting_v28_20210916.rds")) %>% 
  # Keep only US values
  filter(ccode1_COW == "USA") %>% 
  # rename variables
  rename(UNGA_affin_US = UNGA_affinity) %>% 
  # Select variables
  select(ccode2_COW, year, UNGA_affin_US) %>% 
  filter(!is.na(ccode2_COW)) %>% 
  filter(ccode2_COW != "") %>% 
  # Filter years
  filter(year >= 1980)
```


Fifth, information on UN Security Council Membership.

```{r}
# UNSC
dt_unsc <- read.csv(paste0(dir_data, "/DreherSturmVreeland2009_UNSCMembership.csv"),
                        header = TRUE, fileEncoding = "UTF-8-BOM")
```


Finally, data on forest policies from Wuepper et al. 2024. (https://zenodo.org/records/10842614)


```{r}
dt_Wuepper_raw <- haven::read_dta(paste0(dir_data, 
                               "/Wuepper2024_POLICY DATABASE 1.0.dta")) %>% 
  # fix year mistakes
  mutate(Initiated = case_when(Initiated == 2996 ~ 2016,
                               Initiated == 2914 ~ 2014,
                               Initiated == 2104 ~ 2014,
                               TRUE ~ Initiated)) %>%
  # fix country mistakes
  mutate(ISO3 = if_else(ISO3 == "SAF", "ZAF", ISO3)) %>%
  mutate(ISO3 = if_else(ISO3 == "SNV", "SVN", ISO3)) %>% 
  # fix LEX-FAOC143209 (the only law for Kosovo, wrongly attributed to Serbia)
  mutate(ISO3 = if_else(Link == "LEX-FAOC143209", "XKX", ISO3))

dt_Wuepper <- dt_Wuepper_raw %>% 
  filter(Keyword1 %in% c("Biodiversity", "Forest", "Landuse")) %>% 
  group_by(ISO3, Initiated) %>% 
  summarize(n_forpol = n(), .groups = "keep") %>% 
  ungroup() %>% 
  rename(Country.Code_WDI = ISO3,
         year = Initiated)
```



### Merge

```{r}
dt_reg <- dt_sample_1pct %>% 
  # Attach DV, lead by 1 period
  left_join(dt_tcloss_WRI, by = c("Country.Name" = "Country.Name_WDI",
                                  "year" = "L1.year")) %>%
  
  # Attach IMF conditionality
  left_join(dt_INET_raw, by = c("Country.Code" = "ccode_INET", 
                                "year")) %>%
  mutate(IMFBA1 = replace_na(IMFBA1, 0),
         IMFBA2 = replace_na(IMFBA2, 0)) %>%  
  # Attach IMF program dummy
  left_join(dt_IMF55_raw, by = c("cname.IMF55" = "cname",
                                 "year")) %>%
  mutate(IMF55 = replace_na(IMF55, 0)) %>%
  
  # Attach controls
  # WEO
  left_join(dt_WEO, by = c("Country.WEO",
                                 "year" = "year")) %>%
  
  # WDI
  left_join(dt_WDI_var, by = c("Country.Code", "Country.Name",
                              "year")) %>%
  
  # VDem
  left_join(dt_vdem, by = c("ccode_VDem" = "country_text_id",
                                 "year" = "year")) %>% 

  # UNGA
  left_join(dt_unga, by = c("COW.Code" = "ccode2_COW",
                            "year" = "year")) %>%
  # UNSC
  left_join(dt_unsc, by = c("Country.Code" = "ccode",
                            "year" = "year")) %>%
  rename(UNSC_member = unsc) %>% 
  mutate(UNSC_member = ifelse(is.na(UNSC_member), 0, UNSC_member)) %>%
  
  # Attach Wuepper
  left_join(dt_Wuepper, by = c("Country.Code" = "Country.Code_WDI",
                                   "year" = "year")) %>% 
  mutate(n_forpol = replace_na(n_forpol, 0))
```


```{r}
saveRDS(dt_reg, paste0(dir_data, "/dt_IMFFOR_Dec2025_final.rds"))
```


# Regression analysis

Let's load the data.

```{r}
dt_reg <- readRDS(paste0(dir_data, "/dt_IMFFOR_Dec2025_final.rds")) %>%
  # Define factors
  mutate(ccode_FE = as.factor(Country.Code),
         year_FE = as.factor(year)) %>%
  # Subset to 2000-2020
  filter(year >= 2000 & year <= 2019) %>% 
  
  # Impute political stability for 2001 (average 2000 and 2002)
  group_by(Country.Code) %>%
  mutate(polstab_WGI_imp = ifelse(is.na(polstab_WGI) & year == 2001,
                 mean(polstab_WGI[year %in% c(2000, 2002)], na.rm = TRUE),
                 polstab_WGI)) %>% 
  ungroup() %>%
  
  # Drop if missing ccode_WDI
  filter(!is.na(Country.Code))
```


```{r}
# Full observations
dt_sample <- dt_reg %>% 
  filter(!is.na(umd_lnnettcloss_ha)) %>% 
  filter(!is.na(IMF55)) %>% 
  filter(!is.na(lngdp_WDI)) %>% 
  filter(!is.na(gdpgrowth_WDI)) %>% 
  filter(!is.na(cabal_WEO)) %>% 
  filter(!is.na(debt_export_WDI)) %>% 
  filter(!is.na(popgrowth_WDI)) %>% 
  filter(!is.na(forestrent_WDI)) %>%   
  filter(!is.na(agr_GDP_WDI)) %>% 
  filter(!is.na(n_forpol)) %>% 
  filter(!is.na(polstab_WGI_imp)) %>% 
  filter(!is.na(libdem_VDem))
```



## Summary statistics:

```{r}
list_covariates <- c("umd_lnnettcloss_ha", "IMF55", "BA2ENV", "BA2FISCAL", "BA2EXT", "BA2FIN", "BA2NEW", "IMFBA2", "lngdp_WDI", "gdpgrowth_WDI", "cabal_WEO", "debt_export_WDI", "popgrowth_WDI", "forestrent_WDI", "agr_GDP_WDI", "n_forpol", "polstab_WGI", "libdem_VDem", "infl_WDI", "elec_VDem", "UNSC_member", "UNGA_affin_US")

# select variables
dt_table <- dt_reg %>%
  filter(Country.Name %in% unique(dt_sample$Country.Name)) %>% 
  filter(!is.na(umd_lnnettcloss_ha)) %>% 
  filter(!is.na(IMF55)) %>% 
  select(all_of(list_covariates))

# use modelsummary
modelsummary::datasummary_skim(dt_table, fmt = 2, histogram = FALSE,
                               output = "markdown")
```



```{r}
# How many countries in the sample?
# 109 with data
n_distinct(dt_reg$Country.Name)
n_distinct(dt_sample$Country.Name)
```


## Baseline analysis

```{r}
# twfe full model
m_twfe <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
m_infl <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
m_elec <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
m_unsc <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
m_unga <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```


```{r}
models <- list(
  "Forest loss t+1 (ln)" = m_twfe,
  "Forest loss t+1 (ln)" = m_infl,
  "Forest loss t+1 (ln)" = m_elec,
  "Forest loss t+1 (ln)" = m_unsc,
  "Forest loss t+1 (ln)" = m_unga
)

coefs_model <- c("IMF55" = "IMF program", # IMFBA1 IMFBA2 IMF55
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "Table2_Dec25.docx")
```


## Calculate effect

Baseline model: 9.2%

```{r}
exp(0.088) - 1
```

What do our findings in Table 2 mean in substantive terms? Assuming an average three-year duration of an IMF program, our estimates imply that each IMF loan is, at the average, associated with a decrease in forest area of 396km2, almost the size of the Barbados. For illustrative purposes, in our sample of 109 countries between 2000 and 2020, we record 604 country-years with IMF programs. According to our statistical estimates, the marginal increase in tree cover loss associated with IMF programs in these country-years totals 41,130km2, or 2,050km2 per year—about 6% of the global tropical primary forest loss of 2023 (36).


```{r}
summary(dt_sample$umd_tcloss_ha)
```

Mean tree cover loss, in ha: 143,498 = 1434,98 sqkm

1434.98 * 0.092 = 132.0182 per year

396.0545 = per three-year program


```{r}
dt_sample %>% 
  filter(IMF55 == 1)
```

604 country-years with an IMF program

Now total marginal loss:

```{r}
# keep only IMF program years
imf_years <- dt_sample[dt_sample$IMF55 == 1, ]

# marginal effect per country-year (in hectares)
imf_years$me_ha <- 0.092 * imf_years$umd_tcloss_ha

# total marginal forest loss across all IMF country-years (ha)
total_me_ha <- sum(imf_years$me_ha, na.rm = TRUE)

# optional: convert to km^2 (1 km^2 = 100 ha)
total_me_km2 <- total_me_ha / 100
```

```{r}
total_me_km2
```


### F3. IMF programs and deforestation: probing mechanisms

Figure 3. IMF programs and deforestation: probing mechanisms; code developed with ChatGPT.

[first need to run the models in Appendix E]

```{r}
# 1) Put models into named lists by panel
mods_A <- list( # Environmental
  "M1"  = env1,
  "M2" = env2,
  "M3" = env3,
  "M4" = env4,
  "M5" = env5
)

mods_B <- list( # Fiscal issues
  "M1"  = fp1,
  "M2" = fp2,
  "M3" = fp3,
  "M4" = fp4,
  "M5" = fp5
)

mods_C <- list( # External sector
  "M1"  = ext1,
  "M2" = ext2,
  "M3" = ext3,
  "M4" = ext4,
  "M5" = ext5
)

mods_D <- list( # Monetary policy
  "M1"  = fin1,
  "M2" = fin2,
  "M3" = fin3,
  "M4" = fin4,
  "M5" = fin5
)
```



Extract the relevant information.

```{r}
# 2) Helper to extract selected terms
tidy_terms <- function(model_list, terms_keep, panel_title, policy_pretty = NA_character_) {
  map2_dfr(model_list, names(model_list), \(m, lab) {
    broom::tidy(m, conf.int = TRUE, vcov = ~ ccode_FE) |>
      filter(term %in% terms_keep) |>
      mutate(model = lab)
  }) |>
    mutate(
      panel       = panel_title,
      term_pretty = dplyr::recode(term,
        IMF55 = "IMF program",
        .default = ifelse(!is.na(policy_pretty), policy_pretty, term)
      )
    ) |>
    select(panel, model, term, term_pretty, estimate, std.error, conf.low, conf.high)
}
```


Build a panel using the function above.

```{r}
# 3) Build panel-wise data frames
df_A <- tidy_terms(mods_A, terms_keep = c("IMF55", "BA2ENV"),
                   panel_title = "A. Environmental reforms",
                   policy_pretty = "Environmental reforms")

df_B <- tidy_terms(mods_B, terms_keep = c("IMF55", "BA2FISCAL"),
                   panel_title = "B. Fiscal issues",
                   policy_pretty = "Fiscal issues")

df_C <- tidy_terms(mods_C, terms_keep = c("IMF55", "BA2EXT"),
                   panel_title = "C. External sector",
                   policy_pretty = "External sector")

df_D <- tidy_terms(mods_D, terms_keep = c("IMF55", "BA2FIN"),
                   panel_title = "D. Monetary policy",
                   policy_pretty = "Monetary policy")
```


```{r}
# 4) Bind all for plotting
df_all <- bind_rows(df_A, df_B, df_C, df_D)

# Ensure consistent model ordering everywhere
model_levels <- names(mods_A)
df_all <- df_all |> mutate(model = factor(model, levels = model_levels))
```


Now let's visualize this.

```{r}
# ==== FINAL COMBINED FIGURE ====

# ---- 1. Prepare factor levels ----
pretty_labels <- c(
  IMF55    = "IMF program",
  BA2ENV   = "Environmental reforms",
  BA2FISCAL   = "Fiscal issues",
  BA2EXT   = "External sector",
  BA2FIN  = "Monetary policy"
)

df_all_plot <- df_all %>%
  mutate(
    term_pretty = recode(term, !!!pretty_labels),
    term_pretty = factor(term_pretty,
                         levels = c("IMF program",
                                    "Environmental reforms",
                                    "Fiscal issues",
                                    "External sector",
                                    "Monetary policy"))
  )

# ---- 2. Global y-limits ----
global_range <- range(df_all_plot$conf.low, df_all_plot$conf.high, na.rm = TRUE)
pad <- diff(global_range) * 0.1
y_limits <- c(global_range[1] - pad, global_range[2] + pad)

# ---- 3. Function for each panel ----
make_panel <- function(df_panel, title, ylab = "", show_legend = FALSE) {
  ggplot(df_panel,
         aes(x = model, y = estimate,
             ymin = conf.low, ymax = conf.high,
             shape = term_pretty, color = term_pretty,
             group = term_pretty)) +
    geom_hline(yintercept = 0, linewidth = 0.4, linetype = "dashed") +
    geom_errorbar(width = 0,
                  position = position_dodge(width = 0.5)) +
    geom_point(size = 2.6,
               position = position_dodge(width = 0.5)) +
    coord_cartesian(ylim = y_limits) +
    scale_shape_manual(values = c("IMF program" = 16,   # circle
                                  "Environmental reforms" = 17,
                                  "Fiscal issues" = 17,                                 
                                  "External sector" = 17,
                                  "Monetary policy" = 17)) +
    scale_color_manual(values = c("IMF program" = "gray40",
                                  "Environmental reforms" = "black",
                                  "Fiscal issues" = "black",
                                  "External sector" = "black",
                                  "Monetary policy" = "black")) +
    labs(x = NULL, y = ylab, shape = "Coefficient on:", color = "Coefficient on:",
         title = title) +
    theme_minimal(base_size = 11) +
    theme(
      panel.grid.minor = element_blank(),
      axis.text.x = element_blank(),   # remove x-axis labels
      axis.ticks.x = element_blank(),
      legend.position = if (show_legend) "bottom" else "none",
      plot.title = element_text(face = "bold", size = 10, family = "Arial")
    )
}
```


```{r}
# ---- Function for each panel (per-panel y-limits) ----
make_panel <- function(df_panel, title, ylab = "", show_legend = FALSE, y_pad = 0.10) {

  # compute local y-range from this panel's CIs (+ padding)
  rng <- range(df_panel$conf.low, df_panel$conf.high, na.rm = TRUE)
  pad <- diff(rng) * y_pad
  y_limits <- c(rng[1] - pad, rng[2] + pad)

  ggplot(df_panel,
         aes(x = model, y = estimate,
             ymin = conf.low, ymax = conf.high,
             shape = term_pretty, color = term_pretty,
             group = term_pretty)) +
    geom_hline(yintercept = 0, linewidth = 0.2, linetype = "dashed") +
    geom_errorbar(width = 0, position = position_dodge(width = 0.4)) +
    geom_point(size = 1.4, position = position_dodge(width = 0.4)) +
    coord_cartesian(ylim = y_limits) +  # <- panel-specific limits
    scale_shape_manual(values = c("IMF program" = 16,
                                  "Environmental reforms" = 17,
                                  "External sector" = 17,
                                  "Fiscal issues" = 17,
                                  "Monetary policy" = 17)) +
    scale_color_manual(values = c("IMF program" = "gray40",
                                  "Environmental reforms" = "black",
                                  "External sector" = "black",
                                  "Fiscal issues" = "black",
                                  "Monetary policy" = "black")) +
    labs(x = NULL, y = ylab, shape = "Coefficient on:", color = "Coefficient on:",
         title = title) +
    theme_minimal(base_size = 10) +
    theme(
      panel.grid.minor = element_blank(),
      axis.text.x = element_text(size = 6, family = "Arial"),
      axis.ticks.x = element_blank(),
      legend.position = if (show_legend) "bottom" else "none",
      plot.title = element_text(face = "bold", size = 10, family = "Arial")
    )
}
```


```{r}
# ---- 4. Build each panel ----
pA <- make_panel(filter(df_all_plot, panel == "A. Environmental reforms"),
                 "A. Environmental reforms", ylab = "")
pB <- make_panel(filter(df_all_plot, panel == "B. Fiscal issues"),
                 "B. Fiscal issues", ylab = "")
pC <- make_panel(filter(df_all_plot, panel == "C. External sector"),
                 "C. External sector", ylab =)
pD <- make_panel(filter(df_all_plot, panel == "D. Monetary policy"),
                 "D. Monetary policy", ylab = "")
```



```{r}
model_labels <- c("M1",
                  "M2",
                  "M3",
                  "M4",
                  "M5")

legend_panel_E <- df_all_plot %>%
  filter(panel == "D. Monetary policy") %>%
  mutate(
    # For the legend, label the policy term generically as "Policy area"
    term_legend = ifelse(term_pretty == "Monetary policy", "Policy area", as.character(term_pretty))
  ) %>%
  ggplot(aes(x = model, y = estimate,
                               shape = term_legend, color = term_legend)) +
  # Invisible points to trigger legend while keeping plot blank
  geom_point(aes(y = -999), size = 3) +
  coord_cartesian(ylim = y_limits) +
  scale_x_discrete(labels = model_labels, drop = FALSE) +
  scale_shape_manual(values = c("IMF program" = 16, "Policy area" = 17),
                              name = "Coefficient on:") +
  scale_color_manual(values = c("IMF program" = "gray40", "Policy area" = "black"),
                              name = "Coefficient on:") +
  labs(x = NULL, y = " ") +
  theme_minimal(base_size = 10) +
  theme(
    panel.grid = element_blank(),
    axis.text.y  = element_text(color = "white"),
    axis.ticks.y = element_line(color = "white"),
    axis.title.y = element_text(color = "white"),
    axis.text.x  = element_blank(),
    legend.position = "bottom",
    legend.title = element_text(size = 8, family = "Arial")
  )
```


```{r}
# Assemble rows + perfectly aligned bottom labels + legend
row1 <- plot_grid(pA, pB, ncol = 2, align = "hv")
row2 <- plot_grid(pC, pD, ncol = 2, align = "hv")
row_legend <- plot_grid(legend_panel_E, ncol = 1, align = "hv")
row_legend <- row_legend + theme(plot.margin = margin(t = 20, r = 0, b = 0, l = 0))

combined <- plot_grid(row1, row2, row_legend,
                      ncol = 1, rel_heights = c(1, 1, 0.35), align = "v")

final_plot <- plot_grid(combined, legend, ncol = 1, rel_heights = c(1, 0.12))
```


```{r}
final_plot
```


```{r}
# ---- Save at 1.5 columns wide (114 mm) and high DPI ----
ggsave(
  filename = file.path(dir_data, "/Figure3_Dec25.jpg"),
  plot = final_plot,
  width = 114,          # 1.5 columns
  height = 95,          
  units = "mm",
  dpi = 720,
  device = ragg::agg_jpeg,
  bg = "white"
)
```



## F4. Leave one out

Figure 4. IMF programs and deforestation: leaving out one country at a time

Here, we check whether these results are driven by any individual country.

```{r}
# Create an empty dataframe to store results
results <- data.frame(
  Country_Excluded = character(0),
  IMF55_Coefficient = numeric(0),
  Standard_Error = numeric(0)
)

# List of unique countries
unique_countries <- unique(dt_sample$ccode_FE)

# Loop through each country and estimate the model
for (ccode in unique_countries) {
  # Exclude the current country from the data
  dt_sub <- dt_sample[dt_sample$ccode_FE != ccode, ]

  # Estimate the model
  model <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sub, cluster = "ccode_FE")

  # Extract coefficient estimate, standard error
  coef_summary <- summary(model)
  coef_IMF55 <- coef_summary$coefficients["IMF55"]
  se_IMF55 <- coef_summary$se["IMF55"]

  # Store the results in the dataframe
  results <- rbind(results, data.frame(
    Country_Excluded = ccode,
    IMF55_Coefficient = coef_IMF55,
    Standard_Error = se_IMF55
  ))
}
```


```{r}
# Add difference from full-sample coefficient
results <- results %>%
  mutate(deviation = abs(IMF55_Coefficient - 0.08804705))

# Tag most influential observation: KGZ
results <- results %>%
  arrange(desc(deviation)) %>%
  mutate(label = ifelse(Country_Excluded %in% c("KGZ"), as.character(Country_Excluded), NA))
```


```{r}
m_twfe$coefficients["IMF55"]
# 0.08804705
```

```{r}
# Print and save the results dataframe
p4 <- results %>% 
  ggplot(aes(x = IMF55_Coefficient)) +
  geom_density(fill = "skyblue") +
  geom_vline(xintercept = 0.08804705, color = "black", linewidth = 0.75, linetype = "dashed") +
  geom_text(aes(label = label, y = 12), na.rm = TRUE, size = 2) +
  hrbrthemes::theme_ipsum_rc(base_family = "Arial") +
  theme(legend.position = "None",
        axis.title.x = element_text(size = 9),
        axis.title.y = element_text(size = 9),
        axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 8),
        # CRITICAL: shrink outer plot margins
        plot.margin = margin(t = 4, r = 4, b = 2, l = 6)) +
  labs(x = "Coefficient IMF program",
       y = "Density") +
  scale_x_continuous(minor_breaks = seq(0.07, 0.1, by = 0.005),
                     breaks = c(0.07, 0.08, 0.09, 0.10),
                     limits = c(0.070, 0.100))

# Save at journal size: 1.5 columns = 114 mm wide
ggsave(
  filename = file.path(dir_data, "/Figure4_Dec25.jpg"),
  plot = p4,
  width = 114,
  height = 71.25,
  units = "mm",
  dpi = 720,
  device = ragg::agg_jpeg,
  bg = "white"
)
```



```{r}
results %>% 
  arrange(desc(IMF55_Coefficient)) %>% 
  head()
```

```{r}
results %>% 
  arrange(desc(IMF55_Coefficient)) %>% 
  tail()
```


```{r}
# Define selected countries
top_countries <- c("SEN", "EGY", "PER", "KGZ", "TJK", "SLE")

# Summarize by country
country_summary <- dt_sample %>%
  filter(ccode_FE %in% top_countries) %>%
  group_by(ccode_FE) %>%
  summarise(
    IMF_participation_rate = mean(IMF55, na.rm = TRUE),
    Avg_treecover_loss = mean(umd_tcloss_ha, na.rm = TRUE),
    Total_treecover_loss = sum(umd_tcloss_ha, na.rm = TRUE),
    N_years = n()
  ) %>%
  arrange(desc(IMF_participation_rate))
```


```{r}
country_summary
```




# Appendix

## C. Variation in tree cover loss

### Figure C1.

In Figure C1, we visualize the level of annual tree cover loss (non-related to fires) in our sample. Each line represents a low- or middle-income country. As is evident, our dependent variable varies considerably more between countries than within countries over time.

Note that `year` refers to `umd_lnnettcloss_ha` at  t+1, i.e., the data covers the span from 2001 until 2020 (as a function of IMF programs from 2000 to 2019).


```{r}
# Visualize within-country variation
dt_sample %>%
  ggplot(aes(x = year + 1, y = umd_lnnettcloss_ha, color = ccode_FE)) +
  geom_line() +
  hrbrthemes::theme_ipsum_rc(base_family = "Arial") +
  theme(legend.position = "None",
        axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14)) +
  labs(y = "Forest loss ha (ln)",
       x = "Year") + 
  scale_y_continuous(breaks = c(0, 5, 10, 15),
                     limits = c(0, 15.2)) +
  scale_x_continuous(breaks = c(2001, 2010, 2020),
                     limits = c(2001, 2020))
ggsave(filename = paste0(dir_data, "/FigureC1_Variation_Dec25.jpg"),
       dpi = 300, device = "png", height = 5, width = 8, bg = "white")
```



## D. Additional analyses

### D1. Predicting tree cover loss without IMF treatment

Table D1. Predicting tree cover loss without IMF treatment

```{r}
# twfe full model
m_twfe <- feols(umd_lnnettcloss_ha ~ lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
m_infl <- feols(umd_lnnettcloss_ha ~ lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
m_elec <- feols(umd_lnnettcloss_ha ~ lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
m_unsc <- feols(umd_lnnettcloss_ha ~ lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
m_unga <- feols(umd_lnnettcloss_ha ~ lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```


```{r}
models <- list(
  "Forest loss t+1 (ln)" = m_twfe,
  "Forest loss t+1 (ln)" = m_infl,
  "Forest loss t+1 (ln)" = m_elec,
  "Forest loss t+1 (ln)" = m_unsc,
  "Forest loss t+1 (ln)" = m_unga
)

coefs_model <- c("lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableD1_Dec25.docx")
```



### D2. IMF programs (binding reforms) and forest cover

Table D2. IMF programs (binding reforms) and forest cover

```{r}
# twfe full model
m_twfe <- feols(umd_lnnettcloss_ha ~ IMFBA2 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
m_infl <- feols(umd_lnnettcloss_ha ~ IMFBA2 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
m_elec <- feols(umd_lnnettcloss_ha ~ IMFBA2 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
m_unsc <- feols(umd_lnnettcloss_ha ~ IMFBA2 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
m_unga <- feols(umd_lnnettcloss_ha ~ IMFBA2 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```


```{r}
models <- list(
  "Forest loss t+1 (ln)" = m_twfe,
  "Forest loss t+1 (ln)" = m_infl,
  "Forest loss t+1 (ln)" = m_elec,
  "Forest loss t+1 (ln)" = m_unsc,
  "Forest loss t+1 (ln)" = m_unga
)

coefs_model <- c("IMFBA2" = "IMF program (binding condition)",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableD2_Dec25.docx")
```


### D3. IMF programs and forest cover: Different samples

Table D3. IMF programs and forest cover: Different samples


```{r}
# Define countries to be included in the sample
# Forest cover as a share of land area in 2000
dt_sample_2000 <- dt_reg %>% filter(year == 2000)

vec_5pct <- dt_sample_2000$Country.Name[dt_sample_2000$forestarea_shr_WDI >= 5]
vec_10pct <- dt_sample_2000$Country.Name[dt_sample_2000$forestarea_shr_WDI >= 10]
vec_15pct <- dt_sample_2000$Country.Name[dt_sample_2000$forestarea_shr_WDI >= 15]
vec_30pct <- dt_sample_2000$Country.Name[dt_sample_2000$forestarea_shr_WDI >= 30]
vec_50pct <- dt_sample_2000$Country.Name[dt_sample_2000$forestarea_shr_WDI >= 50]

dt_sample_5pct <- dt_sample %>% filter(Country.Name %in% vec_5pct)
dt_sample_10pct <- dt_sample %>% filter(Country.Name %in% vec_10pct)
dt_sample_15pct <- dt_sample %>% filter(Country.Name %in% vec_15pct)
dt_sample_30pct <- dt_sample %>% filter(Country.Name %in% vec_30pct)
dt_sample_50pct <- dt_sample %>% filter(Country.Name %in% vec_50pct)
```


```{r}
# Forest share (2000) thresholds
m_5pct <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_5pct, cluster = "ccode_FE")

m_10pct <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_10pct, cluster = "ccode_FE")

m_15pct <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_15pct, cluster = "ccode_FE")

m_30pct <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_30pct, cluster = "ccode_FE")

m_50pct <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_50pct, cluster = "ccode_FE")
```


```{r}
models <- list(
  "Forest loss t+1 (ln) 5%" = m_5pct,
  "Forest loss t+1 (ln) 10%" = m_10pct,
  "Forest loss t+1 (ln) 15%" = m_15pct,
  "Forest loss t+1 (ln) 30%" = m_30pct,
  "Forest loss t+1 (ln) 50%" = m_50pct
)

coefs_model <- c("IMF55" = "IMF program",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableD3_Dec25.docx")
```



### D4. IMF programs and forest cover: subset by income group

Table D4. IMF programs and forest cover: subset by income group

```{r}
# Define subsets by income group
dt_sample_low <- dt_sample %>% filter(Income.2000 == "L")
dt_sample_lm <- dt_sample %>% filter(Income.2000 == "LM")
dt_sample_um <- dt_sample %>% filter(Income.2000 == "UM")
```


```{r}
# Income levels (as per 2000 WB classification)
m_low <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_low, cluster = "ccode_FE")

m_lowmid <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_lm, cluster = "ccode_FE")

m_upmid <- feols(umd_lnnettcloss_ha ~ IMF55 + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample_um, cluster = "ccode_FE")
```


```{r}
models <- list(
  "Forest loss t+1 (ln) L" = m_low,
  "Forest loss t+1 (ln) LM" = m_lowmid,
  "Forest loss t+1 (ln) UM" = m_upmid
)

coefs_model <- c("IMF55" = "IMF program",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableD4_Dec25.docx")
```


## E. Mechanisms

### ENV

```{r}
# twfe full model
env1 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2ENV + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
env2 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2ENV + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
env3 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2ENV + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
env4 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2ENV + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
env5 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2ENV + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```

```{r}
models <- list(
  "Forest loss t+1 (ln)" = env1,
  "Forest loss t+1 (ln)" = env2,
  "Forest loss t+1 (ln)" = env3,
  "Forest loss t+1 (ln)" = env4,
  "Forest loss t+1 (ln)" = env5
)

coefs_model <- c("IMF55" = "IMF program",
                 "BA2ENV" = "Environmental reforms",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableE1_Dec25.docx")
```



### FP

```{r}
# twfe full model
fp1 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FISCAL + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
fp2 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FISCAL + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
fp3 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FISCAL + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
fp4 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FISCAL + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
fp5 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FISCAL + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```

```{r}
models <- list(
  "Forest loss t+1 (ln)" = fp1,
  "Forest loss t+1 (ln)" = fp2,
  "Forest loss t+1 (ln)" = fp3,
  "Forest loss t+1 (ln)" = fp4,
  "Forest loss t+1 (ln)" = fp5
)

coefs_model <- c("IMF55" = "IMF program",
                 "BA2FISCAL" = "Fiscal policy reforms",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableE2_Dec25.docx")
```



### EXT

```{r}
# twfe full model
ext1 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2EXT + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
ext2 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2EXT + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
ext3 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2EXT + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
ext4 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2EXT + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
ext5 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2EXT + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```

```{r}
models <- list(
  "Forest loss t+1 (ln)" = ext1,
  "Forest loss t+1 (ln)" = ext2,
  "Forest loss t+1 (ln)" = ext3,
  "Forest loss t+1 (ln)" = ext4,
  "Forest loss t+1 (ln)" = ext5
)

coefs_model <- c("IMF55" = "IMF program",
                 "BA2EXT" = "External sector reforms",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableE3_Dec25.docx")
```



### FIN

```{r}
# twfe full model
fin1 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FIN + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
fin2 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FIN + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
fin3 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FIN + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
fin4 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FIN + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
fin5 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2FIN + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```

```{r}
models <- list(
  "Forest loss t+1 (ln)" = fin1,
  "Forest loss t+1 (ln)" = fin2,
  "Forest loss t+1 (ln)" = fin3,
  "Forest loss t+1 (ln)" = fin4,
  "Forest loss t+1 (ln)" = fin5
)

coefs_model <- c("IMF55" = "IMF program",
                 "BA2FIN" = "Monetary policy reforms",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableE4_Dec25.docx")
```



### NEW (non-core)

```{r}
# twfe full model
new1 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2NEW + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# infl_WDI
new2 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2NEW + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + infl_WDI | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# elec_VDem
new3 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2NEW + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + elec_VDem | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNSC_member
new4 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2NEW + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNSC_member | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")

# UNGA_affin_US
new5 <- feols(umd_lnnettcloss_ha ~ IMF55 + BA2NEW + lngdp_WDI + gdpgrowth_WDI + cabal_WEO + debt_export_WDI + popgrowth_WDI + forestrent_WDI + agr_GDP_WDI + n_forpol + polstab_WGI_imp + libdem_VDem + UNGA_affin_US | 
                 ccode_FE + year_FE,
         data = dt_sample, cluster = "ccode_FE")
```

```{r}
models <- list(
  "Forest loss t+1 (ln)" = new1,
  "Forest loss t+1 (ln)" = new2,
  "Forest loss t+1 (ln)" = new3,
  "Forest loss t+1 (ln)" = new4,
  "Forest loss t+1 (ln)" = new5
)

coefs_model <- c("IMF55" = "IMF program",
                 "BA2NEW" = "Non-core policy reforms",
                 "lngdp_WDI" = "GDP (log)",
                 "gdpgrowth_WDI" = "GDP growth",
                 "cabal_WEO" = "Current account",
                 "debt_export_WDI" = "Gov. debt service",
                 "popgrowth_WDI" = "Population growth",
                 "forestrent_WDI" = "Forest rents",
                 "agr_GDP_WDI" = "Agriculture value added",
                 "n_forpol" = "Forest policies",
                 "polstab_WGI_imp" = "Political stability",
                 "libdem_VDem" = "Liberal democracy",
                 "infl_WDI" = "Inflation",
                 "elec_VDem" = "Election",
                 "UNSC_member" = "UNSC member",
                 "UNGA_affin_US" = "UNGA voting affinity US")

modelsummary(models,
             stars = c('*' = .1, '**' = .05, '***' = .01),
             coef_map = coefs_model,
             output = "TableE5_Dec25.docx")
```

